#ifndef __S_TRIANGLE_NURBS__
#define __S_TRIANGLE_NURBS__

#include "element.h"
#include "edge.h"

class STriangle : public ElementImplementation<STriangle> {
public:
	STriangle(Node* v0, Node* v1, Node* v2, Attribute& attribute, size_t id);
	
	virtual const Node* GetNode(size_t i) const {return vertex[i];}
	virtual Node* GetNode(size_t i) {return vertex[i];}
	virtual size_t SizeNodes() const {return 3;}
	
	template<typename T> void EvaluateImplementation(T& x, const Vector& P, size_t attributeHandle) const;
	template<typename T> void IntegrateImplementation(T& x, const Node* node, size_t attributeHandle) const;

	/////////////////////
	Node* vertex[3];
	SEdge* edge[3];
	STriangle* neighbor;

	bool Init(const std::string& attribute);

	void ComputeBaryCentric(const Vec3n& P, nfloat& a, nfloat& b, nfloat& c, const std::string& attribute) const;
	Mat4x4n Transform(const std::string& attribute) const;
	Mat4x4n Transform2(const std::string& attribute) const;

	Vec3n Normal(const std::string& attribute) const;
	nfloat Area(const std::string& attribute) const;
	Vec3n Center(const std::string& attribute) const;

	virtual Matrix Stiffness(const std::string& attribute) const;
	nfloat Mass(const size_t i, const size_t j, const Matrix& beta, const Mat4x4n& M) const;
	virtual Matrix Mass(const std::string& attribute, const bool lumped) const; 
	
	virtual void Dispatch(Dispatcher* d) {
		d->Visit(this);
	}
private:	
	Matrix DetermineCoefficients3(const Vec2n& p0, const Vec2n& p1, const Vec2n& p2) const;
};

template<typename T>
void STriangle::EvaluateImplementation(T& x, const Vector& P, size_t attributeHandle) const {
	assert(P.size() == 3);

	x = T(0.0);
	
	Vec3n w = math::V3n(P(0), P(1), P(2));
	for(size_t i = 0; i < 3; i++) {
		const T& a = vertex[i]->Get<T>(attributeHandle);
		x += w[i] * a;
	}
}

nfloat f1Temp(nfloat x, nfloat y, nfloat z);

template<typename T>
void STriangle::IntegrateImplementation(T& x, const Node* node, size_t attributeHandle) const {
	const nfloat scale = 0.5;
	static nfloat g[16][3] = {//u, v, w
		{ 0.33333333333333331, 0.33333333333333331, scale * 0.1443156076777862 },
		{ 0.081414823414554124, 0.4592925882927229, scale * 0.095091634267284966 },
		{ 0.4592925882927229, 0.4592925882927229, scale * 0.095091634267284966 },
		{ 0.4592925882927229, 0.081414823414554124, scale * 0.095091634267284966 },
		{ 0.89890554336593786, 0.050547228317031033, scale * 0.032458497623198135 },
		{ 0.050547228317031033, 0.050547228317031033, scale * 0.032458497623198135 },
		{ 0.050547228317031033, 0.89890554336593786, scale * 0.032458497623198135 },
		{ 0.65886138449647969, 0.1705693077517601, scale * 0.1032173705347184 },
		{ 0.1705693077517601, 0.1705693077517601, scale * 0.1032173705347184 },
		{ 0.1705693077517601, 0.65886138449647969, scale * 0.1032173705347184 },
		{ 0.0083947774099572114, 0.72849239295540413, scale * 0.027230314174434864 },
		{ 0.0083947774099572114, 0.26311282963463872, scale * 0.027230314174434864 },
		{ 0.72849239295540413, 0.26311282963463872, scale * 0.027230314174434864 },
		{ 0.72849239295540413, 0.0083947774099572114, scale * 0.027230314174434864 },
		{ 0.26311282963463872, 0.0083947774099572114, scale * 0.027230314174434864 },
		{ 0.26311282963463872, 0.72849239295540413, scale * 0.027230314174434864 }
	};
	
	const Vec3n& a = vertex[0]->Get<Vec3n>("position");
	const Vec3n& b = vertex[1]->Get<Vec3n>("position");
	const Vec3n& c = vertex[2]->Get<Vec3n>("position");
	const nfloat J = 0.5 * length(cross(b - a, c - a));

	const Mat4x4n M1 = Transform2("position");
	
	const Mat4x4n M2 = Transform("position");
	Vec4n at = M2 * math::V4n(a, 1.0);
	Vec4n bt = M2 * math::V4n(b, 1.0);
	Vec4n ct = M2 * math::V4n(c, 1.0);

	Matrix beta = DetermineCoefficients3(math::V2n(at[0], at[1]), math::V2n(bt[0], bt[1]), math::V2n(ct[0], ct[1]));
	Vec3n ai;
	for(size_t k = 0; k < 3; k++) {
		if(vertex[k] == node)
			ai = math::V3n(beta(0, k), beta(1, k), beta(2, k));
	}
	
	x = 0.0;
	for(int i = 0; i < 16; i++) {
		const nfloat& u = g[i][0];
		const nfloat& v = g[i][1];
		const nfloat& w = g[i][2];
		Vec4n p = M2 * (M1 * math::V4n(u, v, 0.0, 1.0));

		const nfloat phi = dot(ai, math::V3n(p[0], p[1], 1.0));

   		x += w * Evaluate<T>(p, attributeHandle) * phi * J;
	}
}

#endif //__S_TRIANGLE_NURBS__
